m 


NPS-MA-92-006 

NAVAL  POSTGRADUATE 

Monterey,  California 


SCHOOL 


AD-A255  187 


IJ  I  iC 

El.yCTE 
SEP  2  5  1992 


A 


f 


DFTS  OK  IRREGULAR  GRIDS: 
THE  ANTERPOLATED  DFT 

by 

Van  Emden  Henson 
Technical  Report  for  Period 
October  1990-March  1992 


Approved  for  public  release;  distribution  unlimited 

Prepared  for:  Navoi  Postgraduate  School 
Monterey,  CA  93943 


92-25841 


NAVAL  POSTGRADUATE  SCHOOL 
MONTEREY,  CA  93943 


Rear  Admiral  R.  W.  West,  Jr.  Harrison  Shull 

Superintendent  Provost 


This  report  was  prepared  in  conjunction  with  research  conducted  for 
the  Naval  Postgraduate  School  and  for  the  Institute  for  Mathematics 
and  Its  Applications.  Funding  was  provided  by  the  Naval 
Postgraduate  School.  Reproduction  of  all  or  part  of  this  report  is 
authorized. 


Prepared  by: 


VAN  EMDEN  HENSON 

Asst.  Professor  of  Mathematics 


Reviewed  by: 


Department  of  Mathematics 


Released  by: 


MARTO 
Dean  df  Research 


UNCLASSIFIED 

SECURITv  ClASSif ‘(.ATiQN  '-CS  PAGE 


REPORT  DOCUMENTATION  PAGE 


RE  S'o. 


f  orrn  Appro\,cd 

ova  ^o  C704  0188 


la  REPORT  SECURII'  CLASSif.CAIlOU 

UNCLASSIFIED 

2a  security  CLASSiEiCATiOM  AUThOR'Tv 

2b  OECLASSIEiCATiOU/ DOVvNGRAD  UG  SCHEDULE 

a  PERFORMING  ORGANIZATION  REPORT  NUMBERiS) 

NPS-MA-92-006 

6a  NAME  OE  PERFORMING  ORGAN. ZATiQN 

6b  C'FF  CE  SVVBO'. 

Naval  Postgraduate  School 

,^^1/  applicable) 

6c  ADDRESS  (Cify.  Srafe.  and  /fPCodp) 

Monterey,  CA  93943 

8a  NAME  OF  F  -'J-JI-.O  SP'jNSO“  N'.S 

Bt>  fif;  c:  S -VH' 

ORGAN12ATON  MSF,  AFOSR, 

(It  aprhiabtp) 

Naval  Postgraduate  School 

Be  ADDRESS  (C'f)  Stale  and /'P  Cede) 

Washington,  D.C.  (NFS) 

Washington,  D.C.  (NFS) 

Monterey,  C.A  9  3  94  3 

i  D-S'3  T/  f)C 

Approved  for  public  release; 
distribution  unlimited 


5  VO'JIIOR  fiG  i.'*  ■  .  •-Rfi'On'  .  V 

XPS-MA-92-006 


Naval  Postgraduate  School 


'D  A  jr, Xify  Srare  a->'J/rC  Ic) 

Mcnterey,  CA  93943 


9  vrr  ’.a'- 

OSNN  Direct  rundinc 


11  titie  (Include  Secutdy  Classification) 

DFTS  on  Iri-Focular  Grids:  The  An. terp.-: 

*  1-  r-  n  ^  ^ 

'2  personal  AjI-OR'Si 

Van  Emden  Henson 

tja  type  oe  report 

'  tb  t  (  .  J  .  RF  D 

M  •.'A-;  (Jf  Rfr- 

'  P-  -rear  \\.i-’b 

Technical  Report 

'POM  10/90  TO  0  3/92 

3  1  M  D  r  c  i 

h  loo: 

1  • 

i  -I 

'6  SUPPLEMEN'ARV  f.OTAlK.iN 

coa-. ) 


SgB  QPfj'. 


^3  S'.  -Gf^  :  '•  <i  r.',-*- ,  r-j  <  :*  .’’-G 

Discrete  Fourier  Transform 
FFT 

Anterpolat ion 


19  ABSiRACf  (Continue  on  tevetse  if  nec.essary  and  identify  by  blotk  no— bn/ 

Abstract.  In  many  instamces  the  discrete  Fourier  transform  {DFT)  is  desired  for  a  data  set 
that  occurs  on  an  irtegoiar  grid.  Commoniy  the  data  are  interpolated  to  a  tegular  gr.d.  and  a  fast 
Fourier  transform  (FFT)  is  then  applied.  drawback  to  this  approach  is  that  typically  the  data  have 
unknown  smoothness  properties,  so  that  the  error  in  the  interpolation  is  unknown. 

An  alternative  method  is  presented,  based  upon  multilevel  integration  techniques  introduced  by 
A.  Brandt.  In  this  approach,  the  kernel,  e"'“‘,  is  interpolated  to  the  irrepuior  grid,  rather  than 
interpolating  the  data  to  the  regular  grid.  This  may  be  accomplished  by  pre-multiplying  the  data  by 
the  adjoint  of  the  interpolation  matrix  (a  process  dubbed  anterpolation),  produang  a  new  regular- 
grid  function,  and  then  applying  a  standard  FFT  to  the  new  function.  Since  the  kernel  is  C~  the 
operation  may  be  carried  out  to  any  preselected  accuracy. 

A  simple  optimization  problem  can  be  solved  to  select  the  problem  parameters  in  an  effiaent  way. 
If  the  requirements  of  accuracy  are  ot  strict,  or  if  a  small  bandwidth  is  of  interest,  the  method  can 
be  used  in  place  of  an  FFT  even  when  the  data  are  regularly  spaced. 
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Abstract.  In  many  instances  the  discrete  Fourier  transform  (DFT)  is  desired  for  a  data  set 
that  occurs  on  an  irregular  grid.  Commonly  the  data  are  interpolated  to  a  regular  grid,  and  a  fast 
Fourier  transform  (FFT)  is  then  applied.  A  drawback  to  this  approach  is  that  typically  the  data  have 
unknown  smoothness  properties,  so  that  the  error  in  the  interpolation  is  unknown. 

An  alternative  method  is  presented,  based  upon  multilevel  integration  technicpies  introduced  by 
A.  Brandt,  In  this  approach,  the  kernel.  is  interpolated  to  the  irreytifir  grid,  rather  than 

interpolating  the  data  to  the  regular  grid.  This  may  be  accomplished  by  pre-multiplying  the  data  by 
the  adjoint  of  the  interpolation  matrix  (a  process  dubbed  anterpolation ).  producing  a  new  regular- 
grid  function,  and  then  applying  a  standard  FFT  to  the  new  function.  Since  the  kernel  is  the 
operation  may  be  carried  out  to  any  preselected  accuracy. 

A  simple  optimization  problem  can  be  solved  to  select  the  problem  parameters  in  an  efficient  way. 
If  the  requirements  of  accuracy  are  not  strict,  or  if  a  small  bandwidth  is  of  interest,  the  inethud  can 
be  used  in  place  of  an  FFT  even  when  the  data  are  regularly  spaced. 


1.  The  formal  DFT  and  the  ADFT.  The  DI  T  is  defined  as  an  nix'raiii ni 

tliat  maps  a  !'‘n»tli-.V  compltw-valtiod  soqnoiire  {j’.j.j-i . J'v-i]  to  aiM.ii,cr 

.\  ci.imploy;- v.aliK'd  st'qm'iKf  -jj-q.i-i . xv-i}  hy  the  rule 

.v-i 

(  1  )  i.,.  =  ■‘"j  '  .  for  A-  =  0.  1 . V  -  ]. 

J=U 

.\s  df‘fiii«'d  in  (  1  i.  tlm  DFT  is  jierformod  on  liata  that  am  pre.'Uim'd  to  ho  xivfu 
on  a  roifnlar  arid,  yyitb  Ciinstant  .si'acinsr  hetwcon  the  tiata  points.  Fnrl In'rinf'm.  the 

transform  values  . i'.v-i}  arc  also  prosiinied  to  lie  on  a  regular  grid  in  the 

froqimncy  dnmain.  In  many  applications,  however,  tlie  data  for  a  proldem  am  not 
spaced  regularly.  It  is  of  some  intere.st.  then,  to  lietennine  how  a  discrete  Fourier 
transform  may  Ite  comjtuted  for  such  a  tiata  set.  To  i>erforni  this  ci.uiiputaiion.  we 
develoj)  and  implement  in  one  dimension  an  algoritlim  btised  on  inultih'Vel  inte<iratiun 
techni<|Ues  otitlined  liy  .Achi  Ilrandt  ([2].  [l]).  The  meiliod  presenttul  liepe  r.an  .also  he 
th’veloped  for  higher-tlimensional  problems.  One  ai>i)lication  of  litis  t»'clinitju>'  ti'  i-  in 
the  reconstruction  of  images  from  jtrojections  (inverting  the  Hadon  transform). 

To  begin,  it  is  necessary  to  decide  what  is  meant  by  a  Discrete  Fourier  Iraiisform 
for  irregularly  spaced  data.  Therefore,  the  conce|)t  of  a  formal  DFT  is  introduced, 
which  is  defined  as  follows: 

fVtiisider  any  set  of  ,V  ordered  ])oints  in  the  inti-rval  [0.  A’ ).  -atisfviim 


0  <  J-o  <  J"!  •  •  •  <  -ry-i  <  A' 

and  suitpose  a  vector- valued  function  (grid  funititm)  niSj)  is  sjn'cified.  I  he  f,nnal 
DFT  is  d.'fined  as  the  M  =  .\/i  -b  M2  -r  1  (juantities 


•>-( 

(2)  .  w-(  =  21—.  -.\/i  <  /  <  M2  . 

;=U  ■' 

where  /  is  an  integer,  and  .\/i  and  .t/y  are  positive  integers  sjiecifying  the  range  of 
fretpiencies  of  interest.  The  formal  DFT  may  be  thought  of  as  an  approximation  to  a 
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selection  (  — A/i  <  I  <  M2)  of  the  Fourier  coefficients  of  u{x).  In  this  view,  u{x)  should 
be  regarded  as  an  A'-periodic  function  known  only  at  the  grid  points  Xj. 

It  is  desired  that  this  sum  be  calculated  to  a  prescribed  accuracy,  say  (||u||i.  where 
|lu||i  is  the  discrete  Li  norm  ||u||i  =  A'ote  that  any  grid  spacing 

is  allowed  for  the  Xj  (in  particular,  the  spacing  needn't  be  constant),  that  there  is  no 
relationship  between  the  integers  and  A',  and  that  there  is  no  requirement 

that  these  integers  have  any  special  value,  (such  as  being  powers  of  2).  Calculating  the 
sum  in  (2)  directly  would  have  a  computational  cost  of  0(M.\)  operations.  Instead 
of  forming  this  sum  directly,  though,  an  approximation  to  it  will  he  cfunputed.  using 
an  F FT  to  accelerate  the  computation. 

The  procedure  begins  with  the  definition  of  an  auxiliary  grid.  flv. •  covering  the 
range  of  values  [0.  A’).  Let  A’,  he  an  integer,  whose  value  will  he  determined  shortly, 
and  let  the  grid  spacing  h  be  defined  by  h  =  A’/A'..  Then  the  auxiliary  grid  consists 
of  the  points  yr  —  (r  -  1)/;.  for  t  =  1.2 . V.. 

Suppose  that  the  value  of  some  function  on  the  regular  erid  llv,  is  to  he 

interpolated  to  the  gridpoints  Xj  by  Lagrangian  interpolation.  W'e  will  identify  the 
interpolation  hy  specifying  the  degree  of  the  polynomial  to  he  usod.  Thus,  p-dcirret' 
Lagrangian  interpolation  is  computed  using  a  pcilyncjinial  of  (h-gree  p  or  h'ss.  For  f'ach 
Xj.  the  1  nearest  neighhors  on  the  grid  mu't  he  located.  l.'U  t  he>e  jioints  hi> 

designated  . These  points  shouhi  h<'  chosen  modulo  A.  so  that 

a  point  near  the  limits  of  the  interval  /).  A’ )  m;iy  havo  neishhcjrs  near  ih*’  other  end 
of  the  interval.  (This  is  justified  hecause,  as  will  he  ^el■n  shortly,  the  fuiiction  to  he 
interpolated  for  the  fortn.'d  DFT  is  .V-periodic. )  For  eto  h  x,.  the  ji  *  1  l.a”rati2ian 
ititer[)oIation  weights  are  com[)uted  hy 


(d) 


1=  n 


(•c,  -  ; -.1 1 

-  'l-J. 


and  the  interpol.it ion  of  the  function  y  to  the  griiipoints  x  is  ‘.ii\i.n 


iji  X  )  ^  y~'  w .  !  /  I  i;i  V.  .  ,  I  . 
t.  =0 

Letting  <7  be  the  vector  of  function  v.altie.'  and  </  he  the  \ector  ol  i  III  erpol.tl  eil 

values  (j{x.).  the  inter[)olation  may  he  written  in  m.ttrix  form 


where  is  the  A  x  A',  ititerpolat ion  matrix  mappiiiii  a  function  oti  to  the 

gridpoints  {xj}.  The  entries  of  this  mtitrix  are 


V  ,.( X  , )  if  t;{j.  II 1  =  III 
0  c'lse  . 


We  are  now  ready  to  compute  an  aj)proximation  to  eipiation  (2).  Ihe  strategy 
will  he  to  interi)olate.  for  each  w'/.  vahie.s  of  the  kernel  Irom  the  aiixili.iry  crid 

11  That  is.  equation  (2)  is  ai)proximated  hy 

.\  -1 

J=^) 


(4) 


where 


p 

^  ■ 

»4=l' 

Let  the  column  vector  of  exponential  values  e"*'"'*'*  be  designateii  T;.  ami  the  vector 
of  interpolated  kernel  values  be  denoted  c/  Then  this  interpolation  can  be 

written 

n  ■ 

Notice,  however,  that  c>  is  to  be  used  as  the  l’^‘  row  of  the  matrix  •zivinsr  the  kf-rnel 
of  the  summation  in  f  1).  To  comiiute  ;  as  a  ran-  vector,  the  .'uijoini  of  the  matrix 
equation  equation  is  needed,  namely 

Let  US  define  the  .U-vector  u  —  ii!-;).  tlie  .\’-vector  it  =  hIt,).  and  the  matrices  eivifig 
the  keriiel.s  of  e(iuations  (2)  ami  (4)  a.-^  IT  and  U'.  re-pecti veiy.  Let  th»'  .U  a  matrix 
whose  /■*’*  row  consists  of  c”'*'!'’.  j[,p  p,,jj,ts  on  (,(<  de>ianai<'d  IT  i  it  is  useful 
to  observe  that  this  matrix  consists  of  .\/  consecutive  rows  of  th*’  standard  l>  FI’  kenml 
for  a  uniform  grid  with  .V.  points).  Then  the  fortnal  DFT  is  a[i[>ro\im.ifed 

u  =  lit? 

H’ir 

=  ^  ■ 

This  tiotation  can  Ite  siiii])lifj(.(J  slightly  by  deiioting  the  vector  ere,, ted  by  m))l'iplyini;' 
il  by  [ly]'^ .  a.s  ft.  .Siiice  the  matrix  [ly]^  is  the  adjoint  of  the  l.aeraiician  interfiolation 
matrix,  the  jirocess  uf  computing  h  =  [1'']^  il  has  been  dubbed  nut)  rj^nliittun.  1  hen 
the  approximation  to  the  fi.irmal  DFT  is 

('))  it  s:  W  it  . 

which  we  call  the  Ant(  rpolatid  Di>cnti  Fonritr  Tnut'form  fADFT/. 

The  ADIT,  as  ;i  matrix  multiplication.  re(|uires  (b.U.V.  t  i '[)eraiions,  li,  gener.d. 
.V.  will  exceed  .V.  so  as  a  matrix-vector  multi()lication.  the  ADFT  has  no  ;id\aniaee 
over  (2).  If.  however.  .V.  is  selected  appropriately,  the  apprerxirnation  ran  be  computed 
(piite  rapidly.  Let  .\f.  =  max  { .U] ,  .\/j}.  Then  if  .V.  is  selected  such  that  .V.  >  2.\/.. 
and  at  the  same  time  is  a  number  for  which  an  F FT  module  exists,  then  tlie  fa.-'t 
Fourier  transform  can  be  applied  to  compute  tlie  DFT  summation 

rT7’{r/}  =  -^  ^  mc-^"'/'-  for  /  =  -:^  +  1.-4^  +  2....^  . 

‘  ’  r=l}  4.  ^ 

Hecalling  that  h  —  it  may  be  .seen  that  the  DFT  sutnmation  thereh.re  yields 

( l/.V.)u((2x/)/.V  ).  .Multiplying  by  .V.  thus  yields  a  set  of  value-,  th.it  includes,  as  a 
suiiset.  all  the  desired  values  of  iti^’i). 

Computing  the  ADFT .  then,  consists  of  two  phases: 

1.  il  is  computed  from  d  by  anterpolatiuri:  u  = 

2.  il  is  computed  from  a  by  a  Fast  Fourier  Transform. 


2.  Operation  count  for  the  ADFT.  The  cost  of  computing  the  ADFT  con¬ 
sists  of  the  cost  of  computing  the  interpolation  weights,  the  cost  of  computing  the 
vector  u  =  and  the  cost  of  the  FFT  on  -V.  points. 

Computing  the  (p-f-  1)A’  interpolation  weights.  w„(j-j),  by  the  ft)rmula  in  (.1)  is 
the  cost  of  the  computing  the  numerator,  since  the  regular  spacing  on  Sly  means 
that  the  denominators  of  w„(jrj)  are  independent  of  j.  To  compute  the  numerators, 

V 

the  product  {Xj  -  computed  for  each  Xj.  requiring  2/;+  1  operations. 

TTl  =0 

Then  the  n'^  interpolation  weight  can  be  ol)tained  liy  dividing  by  tlie  [)roduft  of 

F  ' 

(r„,  —  with  the  precomputed  denominator  n  requiring  2 

msO 

tn  ^  n 

operations  for  each  of  the  p+  1  weights  associated  with  the  point  Xj.  The  calculation 
of  the  weights  thus  retjuires  0(.V(p-|-  1 ))  operations.  It  is  important,  however,  to  note 
that  the  calculation  of  the  weights  is  dependent  oidy  on  the  relationshi[;>  befncen  the 
gridi)oints  {y^.}  and  {jj}.  and  is  indepfuident  of  the  data  set.  u(x,).  This  means  that 
if  a  known  set  of  gridpoints  and  a  standard  auxiliary  grid  are  to  l)e  used 

re])eatedly.  tlie  int('r[)olat ioti  weights  w,.(j-j)  tnay  Ite  precomputed  and  stored,  and 
iiet'dn't  he  included  in  the  cost  of  the  algorithtn.  This  will  be  assuiiied  to  be  tin’  case. 

The  matrix  is  .V.  x  .V  and  the  data  vector  i7  is  .V  x  1.  so  the  computation  of 

u  —  would  lie  0{.\  y.)  if  [lerforttK'd  as  a  ttiatri.x-vecior  multiplication.  There 

is.  however,  a  much  more  eflicient  method,  fhe  index  table  Hij.ti)  can  be  >tored 
alotig  with  the  interpolation  weights.  For  each  Xj  ami  for  each  n.  the  value  of  k( j.  n) 
is  the  index  of  the  n'^‘  interpolatioti  neighbor  that  is  used  to  interpolate  from 
to  the  gridfioint  Fhe  periodic  nature  of  the  k<>rnel  being  interpolated  means  that 
the  interjiolation  is  always  to  a  gridjioint  x.  in  the  center  of  the  set  of  p  interpolation 
lioints  (as  is  well  ktiown.  [a],  the  I..igrangiati  interjiolat ioti  is  better  behaved  when  this 
i.-i  the  c.'ise).  If  J)  is  0(1(1.  then  Xj  always  li('s  lietween  ‘'"'I  t-Li,-  ''•hile  if  p 

is  even,  tlien  is  the  closest  gridpoitit  oti  il\^  to  Xj. 

( 'otnputing  t  he  v(>ctor  h  i.'  l hen  very  ea.-y.  aiid  ma\'  !)(>  done  in  2.Vi  /)-  1  1 1 qe  ration,'. 
according  to  the  algorithm; 

1.  Initialize  )'/(  //.)  =  0  for  all  //.  E  (  1  S  "  ^  .N.  t 

2.  For  J  =  0.  1 . V  -  1 

For  n  :=  0.  1 . p 

) )  —  It  f  !Fij. 1 )  +  (  J- J  It i  -r  J ) . 

Having  com[iuted  the  values  of  u.  consider  now  the  cost  of  the  F7  7'  [lortion  of 
the  ADFT.  This  is  sitnply  the  cost  of  an  .V. -point  FFT.  In  the  muxt  section  criteria 
for  choosing  .V.  will  be  determined.  For  now  the  only  requirement  is  that  an  FF  T  can 
be  computed  on  a  vector  of  length  .V..  As  such.  .V,  tmist  have  fai  tors  for  wIFk  Ii  F' F  T 
modules  are  available.  For  the  purpose  of  an  operation  count,  however,  it  is  easi('st 
to  assume  that  .V.  is  a  power  of  2.  Indeed,  we  shall  see  that  we  havi'  great  flexibility 
in  our  selection  of  .V.,  and  since  powers  of  2  or  4  produce  the  most  efficient  7's, 
this  is  a  good  assutnption.  In  this  case,  the  cost  of  the  F  F  T  portion  of  the  .\DF  T  is 
0(.V.  log^  .V.). 
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The  costs  of  the  ADFT  can  now  be  computed.  In  terms  of  data  storage  it  requires 
four  arrays.  One  is  the  A'-point  vector  containing  the  input  data,  u.  In  addition,  an 
•V.-point  complo.v  vector  is  required  for  the  input  and  output  of  the  FFT.  A.^.^uIlling 
that  the  weights  are  precomputed  and  stored,  two  au.viliary  arrays  are  iiores.sary.  an 
.V  X  (p+  1)  real  (or  double  precision)  array  holding  the  inter[)olatioii  weights  w,  (j-j), 
and  the  ,V  x  (p  +  1)  integer  array  of  indices.  n(j,  n). 

If  the  operations  of  multiplication  and  addition  are  counted  equally,  and  if  the 
weights  and  indices  are  pre-stored.  the  operation  count  is  f'l.V.  log’  A  -  complex  o[)er- 
ations  for  the  FFT  portion  of  the  algorithm,  where  C\  depends  on  the  choice  of  FFT 
algorithm.  The  computation  of  u  entails  '2.\(p+  1 )  operations  that  are  real  or  com])lex 
according  to  whether  il  is  real  or  complex.  Counting  both  phases  of  the  algorithm, 
the  operation  count  of  the  ADFT  is 

Cl  .V.  logj  .V.  +  •i.Vlp-u  1)  . 

This  should  l)e  compared  with  the  o))eralion  count  of  the  formal  DFT.  whiidi  is 
0{  M  \  ).  The  computation  of  the  ADFT  is  more  efficient  provided  M  is  larger  tiian 
2(p4-  1  )  —  ( .V.  log_,  .V.  )/.V.  a  condition  that  will  generally  occur  in  pr.ictice. 

3.  Error  analysis  for  the  Anterpolated  DFT.  One  of  the  atir.Kiive  fe.i- 
tures  of  the  ADFT  is  that  the  interpolation  is  iierformed  on  the  hi  it:!  L  which  has 
known  smoothness  [iroperties,  rather  than  the  data  set.  which  generally  h.i>  unknown 
smoothness  [iroperties.  Sinc<‘  interi>olation  error  de])ends  on  the  sinoothnc'S  of  the 
inferpnl.iied  function,  the  error  committed  by  using  the  ADFI  i^  rel;iti\ely  e,i-y  to 
analyse. 

Consider  the  ('rror  in  //-degree  I.agr;ingi;in  polynomial  iiiierpol.at ion .  when  the 
iiiterpol.at ion  is  from  set  of  p  +  1  gridpoinis  that  are  equtdly  .'p.aced.  lau  thesi’ 

gridj/oiiits  be  designtitei  i  se-sl . St  ■  \  fuiiction  f\x).  whose  vabies  are  kiiown  at 

the.-e  gridpoint.-'.  i,>  to  be  interpolati-d  to  the  point  x  6  •sn-s,-]-  I-'"’  -^  =  s"  -  where 
li  is  the  grid'it.aciiig.  and  t  €  [d./>).  The  a|q»roximation  to  /(.f,,  -r  Ihi  is  t!,e  v.alue  of 
the  Lagrangian  interpobation  ;>olynomi;d  i  vi -*-/// ). 

p 

1=0 

Defining  -  ti  t  -  1 )( t  -  2 ) . . .  ( f  -  pj.  and  G  [so  - s;  j-  error  in  t he  inter pol.it ion 

is  bounded  by 

( G )  i/i  s"u  +  tl, )  -  /'.,( 6,  A  lh)\  <  — - ^ - ,  , 

tp  +  1  )■ 

where  =  maxj^-,,  |/('’+'>(s)|.  See  [sj.  p.ages  2bt-2  70.  for  a  deriv.ition  of  this 

error  term. 

It  is  tiseful  to  bound  this  error  more  precisely.  To  do  this,  we  examine  the  lielmvior 
of  the  factoritd  polynomial  ~y(t).  This  polynomial  htis  been  well-studied,  and  many 
results  can  be  haund  in  various  numerical  analysis  texts,  ([.a].  [>■,  [O].  [lb),  i  hese 
results,  however,  are  developed  for  the  case  that  i  can  be  anywhere  in  In 

present  case  the  interpolation  is  always  to  the  ccnfcr  subinterval.  1  hus  for  p  odd, 
t  G  ^4^].  while  for  p  even,  either  /  G  [5-  f  1]  "f  ^  ^  •  §]■ 


To  shorten  the  discussion,  assume  that  p  is  odd.  This  is  the  most  common  case, 
where  p  =  1  gives  linear  interpolation  and  p  =  3  give  cubic  interpolation.  Similar 
results  can  be  obtained  for  p  even.  Consider  the  following  lemma,  the  proof  of  which 
may  be  found  in  [0]. 

Lemma  1.  If  p  is  a  positive  odd  integer,  then 


ma.v 


iMOl 


< 


{p+ 1)!  2f+' 


This  result  can  be  used  to  find  an  error  bound  for  tlm  AD  FT.  In  this  case  the 
functions  being  interpolated  are 

for  /=  +  . 

from  this  set  of  functions,  the  only  ones  whose  values  are  of  interest  are  those  for 
/  between  —M\  and  Recalling  the  definition  .1/.  =  ma.\  { .l/j .  .!/>}.  the  largest 

absolute  value  among  the  freqiKuicies  of  interest  is  =  ( 'ir )/.V  .  Therehme 


0'^ 


max 


'  i  ilfV 


\l) 


It  I 


Inserting  this  and  the  bound  from  I.emma  1  into  e(|uation  itii  ei'.es 


(7) 


» i 


r  i; 


which  is  t  hen  used  to  obtain 

1'/  —  U  hi  =  '  ^  XTi'  "'I-'’  ji~  U' II  ;  .  ,r  e 

-  *1' ~  'i 

=  (// -■\f . 


Finally,  substituting  li  =  .V/ .V.  and  =  (2rr.\/.)/.V  establishes  the  desired  error 
bound.  The  error  in  the  .\Dl  I  ap[)roximation  to  the  iormal  Dl-  T  tiouiided.  for 
-.1/,  <  /  <  .\r.  by 


where  =  '2~l/.\.  and  the  .\DFT  (.a)  is  computed  using  an  FFT  of  length  .V.. 
Since  the  bound  holds  for  all  desired  values  of  w-;.  it  then  follows  that 


(0) 


6 


where  ||  •  is  defined  as  the  maximum  absolute  value  in  the  vector.  It  is  also  worth 
noting  that  an  error  bound  for  any  desired  frequency  can  be  obtained  by  replacing 
u,\\i  with  u,';  in  the  derivation,  leading  to 

(10)  <  (^-^y  '  -Vjlnli,  . 

This  is  especially  useful  information  for  those  occasions  when  only  the  low  freijuency 
components  are  of  interest,  or  when  the  accuracy  required  of  the  appro.ximation  is 
f,reater  for  the  low  frequencies  than  for  tlie  liigh  frequencies. 

4.  Selection  of  p  and  .V..  The  error  bound  just  derived  is  useful  in  that  it 
provides  a  way  to  select  the  operational  parameters  .V.  ami  p.  Recall  that  the  goal 
is  to  calculate  an  appro.ximation  to  u  to  .some  prescribed  accuracy.  \  u  -  U  ■fi\  <  f|!a!!. 
In  practice  we  will  want  to  make  the  error  small,  so  it  will  be  a.-^sumed  that  <  1. 
Comparison  with  (s^)  gives  the  requirement 


f 


which  may  he  written  as 


(11) 


.V.  >  .\/.T 


for  a  given  formal  IJ I  T .  the  lalues  of  .V.  and  >  .tr<'  ron.'idered  to  he  part  of 
the  prohleiu  s[)ecificiition.  lo  en>ure  thtit  the  specified  tu'ciiritcy  is  obtained,  it  is  oi.ly 
necessary  to  select  integers  .V.  and  />  so  th.it  i  !  1 )  is  satisfied.  .\.at iiraP v.  there  may  he 
man;e  comhination,s  of  [naratiiefer  v.'dues  thtit  achieve  this  go.al.  1  he  par.itneters  'hoiild 
therefore  he  selected  to  fuillill  scjine  other  desirable  property  as  well.  Specilictdl v.  thev 
should  he  selected  also  to  minimize  the  comput.'ilional  effort  of  the  .algririthm. 

lo  see  how  this  may  he  ;iccom()li.shed.  recall  that  the  work  in'.'i,|\-ed  in  comput¬ 
ing  the  .\Dh  I  with  .N .  points  tin  the  tiu.xiliary  grid  Q\  .and  /(-degree  htigriuigian 
interpolation  is 


C(.V.  log.  .V.)  +  (;(.V(/>-|-  l  i)  . 

The  value  of  the  constant  on  the  1))  term  dejieiids  on  whether  the  weights 

and  indices  are  pre-stored.  or  calculated  "on  the  fly".  For  the  amdysis  that  follows, 
we  assume  the  weights  atul  indices  are  pre-stored.  in  \\hich  ctise  the  constant  is  2. 
The  constant  on  the  first  term  depends  on  si’veral  factors.  I  l  ls  gener.allv  have 
a  coitiple.xity  of  ( .\ /r/)  log(  .\ /r/)  for  some  number  7  >  1.  If  the  data  have  certtun 
symmetries,  then  a  specialized  F  FT  may  be  used  for  faster  computation  (['10  [7]. [10’). 
1  he  variety  of  available  F FT  algorithms  pursuades  us  to  leave  the  constant  on  the 
first  term  as  an  unspecified  i)arameter.  ( '1 . 

The  total  work  iit  computing  the  ADFT  can  therefore  he  written  as  a  function  of 
the  two  parameters  .\.  and  p.  For  a  fi.xed  problem  size  (.V  and  and  a  [trescrihed 
error  tolerance  (,  the  work  in  computing  the  ADFT  to  the  required  accuracy  is 


(12) 


lF(.V.,p)  =  Ci-V.  log^  .V.  +  2.V(p  +  1)  . 


and  we  seek  an  optimal  parameters  minimizing  U'(  .V,,p)  over  all  combinations  ( -V..p) 
satisfying  ( 1 1 ).  if  such  a  choice  exists. 

Limiting  cases  may  be  determined  by  examining  nmnst  itfighbor  inter[)nlation 
(p  =  0),  as  well  as  extremely  high  degrees  of  interpolation  {p  ~  tc).  Substituting  the 
limiting  values  of  p  into  (11),  and  noting  that  equality  will  suflice  to  ensure  that  the 
required  accuracy  is  attained,  we  obtain  bounds  for  the  selection  of  .V.,  namely 


-U.T  <  .V.  < 


■l/.rr.V 

( 


for  all  values  of  p  >  0. 

The  existence  and  uniqueness  of  optimal  solutions  are  fairly  ea.sy  to  establish. 
ir(.V..p)  is  continuous  with  respect  to  each  of  its  variables,  and  both  of  the  first 
partial  derivatives  are  everywhere  positive.  This  observation  leads  to 

Lem.ma  2.  L( t  S  be  tht  set  |(.V.,p):  .\.  >  .1/. ;r(.V/t  [f  / 

that  portion  of  the  boundary  of  S  given  by  .V.  =  .1/.  .V/(  | .  lion 

t/(-fo-yij)  €  S.  there  exists  n  point  if.rj]  G  itS  such  that  U’(  f-n}<  »■(  ■f(i-  ,7u  )• 

Proof:  Since  (j-,). G  -S.  the  i)oint  (s-,'Ai)  G  where  t  j  ’  * ‘  . 

Furthermore.  F  <  j\,.  Then  since  the  partial  derivative  i.if  the  work  functiiui  with 
respect  to  .V.  i,'  everywhere  positive,  <  U’l  xo. //n ).  | 

I  he  utility  rif  Lemma  2  is  tent  the  o,)timization  problem  can  be  reu  riitiui  as  a 
problem  in  a  sin.gle  variabhu  Since  for  every  point  in  S  tlu're  is  some  point  ahuie  US 
that  recjuires  less  work,  it  is  only  n(-cessar\  to  seek  a  n.inimum  from  the  point-  of  i/s. 
This  can  be  done  by  parameti'rizing  .V,  and  ji  as  f-'iiciions  of  a  single  variable. 


(Id) 


Then  on  it.S  we  find  that 

(11)  .\ .  =  .U.  rro  aiid  p  — 


Since  0  <  p  <  >:.  the  value  of  b  is  restricted  t  )  the  int<'rv.il  (F.V/f].  Substituting 
these  exjiressions  into  (12).  the  work  eipiation  may  be  rewritten  as  a  function  of  b 
alone 

( Id)  \V(b)  =  C\M.-b  \n^,(M.-b)  +  2.V  log,  . 

and  the  problem  is  to  minimize  (  Fa)  subject  to  the  constraint  1  <  6  <  i  .V,'(  ).  Once  b 
is  determined,  the  necessary  values  of  .V.  .uid  p  r.an  tie  (ditaim'd  from  f  1  1  i.  W’e  m.ay 
now  ('stal)lish 

Theorem  1.  There  exists  a  nnif/ue  value  b  that  tntnu/ii:i  s  /I'>)  .'•uhjevt  tn  1  < 
b  leq{\/e).  Therefore  the  irork  function 


ir(  -Y..  p)  =  C'l  .Y.  log^,  .Y.  -+  2.Y(p  +  1 ) 

has  a  unique  rninirnuni,  subject  to  the  constreu  its 

y. 

8 


and  0  <  p  <  ac 


Proof:  U’(6)  is  continuous  and  differentiable  with  respect  to  b  on  (LA/r].  Dif¬ 
ferentiating  equation  (15)  yields 


=  A'l  ln(  h  ^b)  — 


bilnbj^ 


where 


A';  =  LlH:!  _  ^  oA/.t  .  and  A, 

In  2 


■JA-ln(f) 


for  U  =  0,  then,  b  must  satisfy  6(ln6)‘lrd A.>i)  =  h'^/ h\.  N()\  ■  ir'  is  al.'O  coiit iii- 
uous  and  differenti.ible  on  ( 1.  A'/c].  and  differentiating  yields 

.  ,■  f  \r)b  +  -2\ 

"  ■ 

Since  b  >  1  we  s('e  tliat  ir"(/j)  >  0  fc)r  all  b  £  |1.A'/(1.  mi  any  criiii'.il  p"i;,r  in 
interval  must  correspond  to  a  local  mini:num. 

It  is  ajipareiit  t liat  -x  a^  6  —  1.  L'xamination  of  tin' endpuini  f  =  A  ,  ' 

re\eals  that  since  A'l  >  I\j  >  t.  a. id  <  >C  1.  we  liave  that  if'i  A'  '<  i  >  II.  1  ui'fer. 
ir"(6)  >  0  implies  that  \\''{bj  h.as  exactly  oi.e  >ign  cliang''  in  tin’  interval  ■  i .  .V  .  , 
The  pc.iint  h^j  at  which  this  occurs  is  therefore  .a  global  tninimu’  -  f  ir  ll'ifi.  ami  the 
value  ir(  A’../i).  where 


I  ■'  \ 

A.  =  M.'bfj  and  />  -  log.  I — j 


is  the  iinifjue  global  minimum  for  U’  cui  i>S.  | 

’I'lte  valui’s  of  A',  and  p  obttiined  in  this  m.anner  tire  re.tl  numbers,  l  leue  i-  onlv  a 
limited  I.  urn  her  of  integers  for  which  effn  ient  /-'A’/  s  exist .  and  Lagr.ingiaii  ii.t  erpol.at  ion 
reijuires  p  to  be  an  integer.  Ftirther.  this  entire  discussion  has  lieen  predic.it.  ,1  on  tin’ 
a.ssumption  that  p  is  an  odd  integer,  alihougb  a  siinii.ir  anaiv'ls  min  be  made  f,,r  p 
even.  Once  liie  tlieoretical  valties  of  A',  am!  pare  determined,  they  must  be  modified 
to  allow  computaticin.  There  is  s  mie  flexibility  iii  this,  but  cert.dtilv  selrcting  .\ .  to 
be  tlie  first  integer  larger  tlian  M.~b  for  whicb  an  III  exists,  and  choosing  p  to  bt- 
the  smallest  (tdd  integer  greater  than 


S) 


will  suflice. 

In  order  to  find  the  optimal  vaimm  of  p  and  .V.  it  is  neie>.^ary  to  find  the  v.aiue  i/f 
b  satisfying 


(17)  h(ln/D‘lii(  A.//I  =  -d  , 

W  hile  an  analytic  solution  of  this  equation  caiinot  be  found.  .Xewton's  iteration  may 
be  used.  Table  1  displays  optimal  parameters  .V.  and  p  for  several  combinations  of 
A’.  A/.,  and  c. 


A' 

M. 

e 

N. 

P 

A 

M. 

£ 

N. 

P 

32 

8 

.1 

48.7 

7.7 

128 

32 

.1 

193 

9.9 

32 

32 

.1 

142.6 

15.5 

128 

64 

.1 

325.7 

13.8 

32 

64 

.1 

257.6 

22.2 

128 

128 

.1 

570.7 

19.4 

32 

8 

.01 

52.9 

9.8 

128 

32 

.01 

206.7 

12.1 

32 

32 

.01 

150.1 

19.1 

128 

64 

.01 

344.1 

16.6 

32 

64 

.01 

267.8 

27.1 

128 

128 

.01 

595.6 

23,1 

Table  1.  Optimal  parameters  N,  and  p  computed  for  various  problems. 

5.  An  .AD FT  Example.  To  illustrate  the  ADFT.  consider  the  problem  of  com¬ 
puting  the  formal  DFT  of  the  function  u{x)  =  [(tr  -  j-)/']*.  sampled  on  an  irregular 
grid.  The  irregular  grid  consists  of  .V  =  12''  jeiints  Xj  randondy  spaced  in  the  interval 
tO,2T).  Since  the  extent  of  the  interval  is  2;r.  the  frequencies  _■/  are  just  the  integers 
1.  and  the  formal  DFT  is 


v-i 

ml)  -  Yi  ii{x.]e-‘‘- .  -G4</<G1. 

j=0 

rii<'  sampled  data  are  ,>ho\vn  in  the  top  of  Figure  1.  The  real  part  of  i!^)  is  plotted 
on  the  bofom  of  Figure  1.  The  .ADFT  was  ti.^ed  to  a])])ruximate  tlu'  formal  DFl  . 
'.Mth  va!ua>,  (jf  .V.  =  12'i,  2'G.  512,  and  1021.  Figure  2  displays,  for  each  choice 
of  .\..  the  alisoluti'  \aliie  of  the  error  :'/(/)  —  dl  'F'Oi-  plotted  as  a  function  of  /. 
Linear  interpolation  ip  =  1)  was  used  in  each  case.  Note  that  increasing  the  value 
of  .V.  [modures  a  iioticeable  ch'crease  in  the  error,  and  that  the  error  increases  with 
inrrea'ing  ssavenumber.  as  might  be  inferred  from  { 10).  Figure  3  displays  the  elfect  of 
Using  dilferent  values  of  p  for  fixed  .V..  It  may  be  seen  that  the  error  decreases  rapidly 
a>  ji  i.'  incre.i.o'd.  F(iuation  >9)  predicts  that  the  error  shotild  decrease  at  least  a.s  Last 

a.'  !  }  d('crease>  as  p  or  .V.  are  increased.  Table  2  gives  both  the  infinity  and 

[.~  i.urms  of  the  eriair  \u{l}  —  [lFiij(/i|  for  several  valtit's  of  p  and  each  of  .5.  =  25G 
and  .V.  =  512.  For  .V.  =  25G.  the  error  bound  decreases  by  O.GIG')  each  time  p  is 
iticrea.sed  by  2.  The  experimental  error  is  diminished  by  a  factor  of  a[)proximately 
0.3  as  })  is  increased  from  1  to  3.  and  by  a  factor  of  approximately  0.1  N'.i'h  each 
'ucreediijg  increase,  better  than  the  theory  predicts.  Similarly,  for  .V.  =  512,  the 
theoretical  bound  decreases  by  0.15121  as  p  is  increased  by  2.  while  the  experimental 
decrease  is  approximately  0.11  for  each  increase,  a  slightly  better  result.  .Numerical 
experiments  on  numerous  other  irregularly  sampled  functions,  with  variotis  degrees 
of  Minoijthness,  produced  similar  results.  In  tiiese  experiments  the  .A DFl  behaved 
in  a  similar  fa.'hicjii  as  it  did  for  the  Linction  discussed  above.  There  is  dramatic 
improvement  with  increasing  values  of  .V..  and  p.  .-\s  might  be  expected,  the  error 
diminished  faster  with  smooth  functions  than  tiiscontinuous  functions. 


■B 

D 

WErrorW^ 

llError]], 

— 

■B 

B 

■lijaxuiw 

256 

1 

1 .20GG3 

0.434861 

■ 

512 

D 

0.290.501 

0.109943 

25G 

3 

0.357S89 

0.116080 

■ 

512 

3 

0.029351 

0.00.8315 

25G 

5 

0.1  44478 

0.039136 

■ 

512 

■a 

0.003360 

0.000817 

256 

< 

0.061305 

0.014983 

■ 

512 

7 

0.000419 

9.6638.5e-05 

Table  2.  Errors  of  the  .ADFT  for  various  values  of  ,V.  and  p. 


10 


6.  Some  Open  Questions  about  the  ADFT.  Like  the  continuous  Fourier 

transform,  the  DFT  has  several  important  properties,  such  as  linearity,  the  convolu¬ 
tion  and  correlation  properties,  the  shifting  property,  the  modulation  property,  and 
Parseval’s  relation.  To  what  extent  these  properties  hold  for  the  ADFT  is  an  open 
question.  The  linearity  holds  can  be  established  immediately,  by  noting  that  both  the 
formal  DFT  and  the  ADFT  can  be  written  as  matrix  operations,  so  they  are  linear 
operators.  Certain  symmetry  properties  are  ea.sy  to  establish.  For  example,  applying 
the  ADFT  to  a  real-valued  vector  will  yield  a  conjugate  symmetric  result,  that  is 
u(w)  =  because  the  vector  is  real-valued,  and  because  the  ADFT  is 

computed  by  applying  the  F FT  operator  to  this  vector.  The  DFT,  and  therefore  the 
F FT ,  maps  a  real  vector  to  a  conjugate  symmetric  vector  [4].  .Applying  the  DFT  to 
data  vectors  with  other  symmetries  (even.  odd.  quarter-wave,  etc.)  yields  output  vec¬ 
tors  with  other  types  of  symmetries  [10].  It  is  natural  to  ask  which  of  these  symmetry 
properties  are  inherited  by  the  formal  DFT  or  the  ADFT.  It  seems  reasonable  to 
postulate  that  if  the  irregular  gridpoints  are  symmetrically  disposed  and  the  function 
u(j-j)  is  symmetric  then  the  symmetry  property  of  the  DFT  might  he  inherited  by 
the  formal  DFT  and  the  .\DFT . 

.An  important  (piestion  is;  How  is  tin*  formal  DFT  related  to  th'>  continuous 
Fourier  transform?  That  is.  to  what  ext^uit.  ;iiid  with  what  i  rror.  iloes  the  forina!  DFl 
approximate  the  /■'T'?  .Answering  this  <juestion  may  [irove  to  he  a  leiigthy  process. 
.M;iny  related  ([uestions  will  also  arise.  For  example,  how  does  the  sampling  theorem 
apply  to  an  irregular  grid?  What  freejuencies  can  be  rej)resent»'d  acniralely.  and  what 
constitutes  aliasittg?  Is  there  some  an.alog  to  the  Foisson  summation  theorem?  .Many 
problems  feature  irreguhirly  spaced  data.  S(,i  it  may  lie  assumed  that  the>e  questions 
are  of  some  interest . 
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Fig,  2.  The  absolute  vajue  of  the  error  of  the  ADFT.  plotted  as  a  function  of  wavenumber.  Each 
graph  corresponds  to  a  different  cnoice  of  .V.,  while  linear  interpolation  'p  =  li  is  used  for  ail. 


Fig.  ■'].  The  absolute  vaJue  of  the  error  of  the  ADFT,  plotted  as  a  function  of  wavenumber.  S,  =  old 
for  each  graph,  while  several  values  of  p  are  used. 
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